#include <cassert>
#include <algorithm>
#include <kuai/mol/ct.hpp>


namespace kuai {

	void ConnectTable::makeBuffer(const Index input[], Index buffer[]) {
		if (input != buffer) {
			memcpy(buffer, input, sizeof(Index) * (2 + input[0] + input[1] * 3));
		}
		if (input[0] > 0) {
			ConnectTable table(buffer);
			for (size_t i = 0; i < buffer[0]; ++i) {
				table._degrees[i] = table._offsets[i] = 0;
			}
			Index* bonds = reinterpret_cast<Index*>(table._bonds);
			for (size_t i = 0; i < buffer[1]; ++i) {
				table._degrees[bonds[i*3+1]] += 1;
				table._degrees[bonds[i*3+2]] += 1;
			}
			table._offsets[0] = 0;
			for (size_t i = 1; i < buffer[0]; ++i) {
				table._offsets[i] = table._offsets[i-1] + table._degrees[i-1] * 2;
				table._degrees[i-1] = 0;
			}
			table._degrees[buffer[0]-1] = 0;
			for (size_t i = 0; i < buffer[1]; ++i) {
				Index atom1 = bonds[i*3+1];
				Index atom2 = bonds[i*3+2];

				table._adjs[table._offsets[atom1] + table._degrees[atom1]*2] = atom2;
				table._adjs[table._offsets[atom1] + table._degrees[atom1]*2+1] = i;
				table._adjs[table._offsets[atom2] + table._degrees[atom2]*2] = atom1;
				table._adjs[table._offsets[atom2] + table._degrees[atom2]*2+1] = i;

				table._degrees[atom1] += 1;
				table._degrees[atom2] += 1;
			}
		}
	}
	std::vector<boost::shared_ptr<ConnectTable::SubMolType> > getConnectedParts(ConnectTable& mol) {
		return ConnectionPartVisitor<ConnectTable>::getConnectedParts(mol);
	}
	std::vector<boost::shared_ptr<ConnectTable::SubMolType> > getConnectedParts(ConnectTable::SubMolType& mol) {
		return ConnectionPartVisitor<ConnectTable>::getConnectedParts(mol);
	}

	Index setupAromaticBond(ConnectTable& table, const Index bondrank[], const Integer charge[]) {
		Index nBonds = table.countBonds();
		std::vector<Index> result(table.countBonds(), 0);
		for (Index i = 0; i < nBonds; ++i) {
			result[i] = table.getBond(i).order;
		}
		Index n = setupAromaticBond(table, bondrank, charge, &result[0]);
		for (Index i = 0; i < nBonds; ++i) {
			table.getBond(i).order = result[i];
		}
		return n;
	}

	void getTotalBondOrder(const ConnectTable& ct, const Integer charges[], Index result[]) {
		Index nAtoms = ct.countAtoms();
		memset(result, 0, sizeof(Index)*nAtoms);
		static const Element& ELEMENT_H = Element::getElement(1);
		for (Index i = 0; i < nAtoms; ++i) {
			const Atom& atomI = ct.getAtom(i);
			Index ndeg = ct.degree(atomI);
			for (Index j = 0; j < ndeg; ++j) {
				const Bond& bondJ = ct.getNeighborBond(atomI, j);
				result[i] += bondJ.order;
			}
			if (charges != NULL) {
				if (atomI.element().elecAffinities < ELEMENT_H.elecAffinities) {
					result[i] += std::abs(charges[i]) * Bond::SINGLE_BOND;
				}
				else {
					result[i] -= charges[i] * Bond::SINGLE_BOND;
				}
			}
		} 
	};
	void getImplicitHydrogen(const ConnectTable& ct, const Integer charges[], Index result[]) {
		getTotalBondOrder(ct, charges, result);
		Index nAtoms = ct.countAtoms();
		for (Index i = 0; i < nAtoms; ++i) {
			const Atom& atomI = ct.getAtom(i);
			result[i] = atomI.element().implicitH - result[i] /  Bond::SINGLE_BOND;
			if (result[i] > 8) {
				result[i] = 0;
			}
		}
	};

	namespace {
	
		static std::vector<RingFinder<ConnectTable>::ResultPointer>
			getAromaticRingCandidates(const ConnectTable& ct, const Index bondrank[],const Integer charges[]) 
		{
			Index nAtoms = ct.countAtoms();
			std::vector<Index> implicitH(nAtoms);
			getImplicitHydrogen(ct, charges, &implicitH[0]);

			std::vector<Atom*> atoms; atoms.reserve(nAtoms);
			for (Index i = 0; i < nAtoms; ++i) {
				const Atom& atomI = ct.getAtom(i);
				switch (atomI.type) {
				case 6: {
						Index deg = ct.degree(atomI);
						if (deg >= 2 && deg + implicitH[i] == 3) {
							atoms.push_back(const_cast<Atom*>(&atomI));
						}
					}
					break;

				case 7: {
						Index deg = ct.degree(atomI);
						if (deg >= 2 && deg + implicitH[i] <= 3) {
							atoms.push_back(const_cast<Atom*>(&atomI));
						}
					}
					break;

				case 8:
				case 16: {
						Index deg = ct.degree(atomI);
						if (deg == 2 && implicitH[i] == 0) {
							atoms.push_back(const_cast<Atom*>(&atomI));
						}
					}
					break;
				}
			}

			Index nBonds = ct.countBonds();

			std::vector<Bond*> bonds; bonds.reserve(nBonds);
			for (Index i = 0; i < nBonds; ++i) {
				if (3 <= bondrank[i] && bondrank[i] < 10) {
					const Bond& bondI = ct.getBond(i);
					if (bondI.order < Bond::TRIPLE_BOND) {
						bonds.push_back(const_cast<Bond*>(&bondI));
					}				
				}
			}

			nAtoms = atoms.size();
			nBonds = bonds.size();

			std::vector<RingFinder<ConnectTable>::ResultPointer>  results;
			
			while (true) {
				SubMolecule<ConnectTable> candidates(const_cast<ConnectTable&>(ct), atoms, bonds, true);				
				results = RingFinder<ConnectTable>::getRingBlocks(candidates, 9);
				
				atoms.clear();
				for (std::vector<RingFinder<ConnectTable>::ResultPointer>::iterator
					it = results.begin(); it != results.end(); ++it)
				{
					Index nAtomI = (*it)->countAtoms();
					for (Index j = 0; j < nAtomI; ++j) {
						Atom& atomJ = (*it)->getAtom(j);
						if (atomJ.type == 6) {
							Index deg = (*it)->degree(atomJ);
							Index k;
							for (k = 0; k < deg; ++k) {
								if ((*it)->getNeighborBond(atomJ, k).order > Bond::SINGLE_BOND) {
									break;
								};								
							}
							if (k == deg) {
								continue;
							}
						}
						atoms.push_back(&atomJ);
					}
				}

				if (nAtoms == atoms.size()) {
					break;
				}
				else {
					nAtoms = atoms.size();
				}
			}

			return results;
		}

		Integer countPiElectron(const SubMolecule<ConnectTable>& ring, const ConnectTable& table, const Integer charges[]) { 
			Integer result = 0;
			Integer formalCharge = 0;
			Index nAtoms = ring.countAtoms();
			for (Index i = 0; i < nAtoms; ++i) { 
				const Atom& atomI = ring.getAtom(i);
				if (charges) {
					formalCharge += charges[table.indexOf(atomI)];
				}
				if (isCarbon(atomI)) { 
					++result;
				}
				else if (isNitrogen(atomI)) { 
					Index degree = ring.degree(atomI);
					Index totalOrder = 0;
					for (Index j = 0; j < degree; ++j) {
						const Bond& bond = ring.getNeighborBond(atomI, j);
						totalOrder += bond.order;
					}
					if (totalOrder == Bond::TRIPLE_BOND) { 
						++result;
					}
					else if (totalOrder == Bond::DOUBLE_BOND) { 
						result += 2;
					}
				}
				else if (isOxygen(atomI) || isSulphur(atomI)) { 
					result += 2;
				}
			}
			return result-formalCharge;
		}
	}

		
	std::vector<RingFinder<ConnectTable>::ResultPointer>
		getAromaticBlocks(const ConnectTable& table, const Index bondrank[], const Integer charges[])
	{
		std::vector<RingFinder<ConnectTable>::ResultPointer>
			candidates = getAromaticRingCandidates(table, bondrank, charges);
		std::vector<RingFinder<ConnectTable>::ResultPointer> result;
		for (std::vector<RingFinder<ConnectTable>::ResultPointer>::const_iterator
			it = candidates.begin(); it != candidates.end(); ++it)
		{
			if (countPiElectron(**it, table, charges) % 4 == 2) {
				result.push_back(*it);
			}
			else if ((*it)->countAtoms() < (*it)->countBonds())	// More than one ring
			{
				std::vector<RingFinder<ConnectTable>::ResultPointer>
					rings = RingFinder<ConnectTable>::getSingleRingsFromRingBlock(**it, 9);
				for (std::vector<RingFinder<ConnectTable>::ResultPointer>::const_iterator 
					jt = rings.begin(); jt != rings.end(); ++jt)
				{
					if (countPiElectron(**it, table, charges) % 4 == 2) {
						result.push_back(*jt);
					}
				}
			}
		}
		return result;
	}

	Index setupAromaticBond(const ConnectTable& table, const Index bondrank[], const Integer charges[], Index result[]) {
		std::vector<RingFinder<ConnectTable>::ResultPointer>
			blocks = getAromaticBlocks(table, bondrank, charges);
		for (std::vector<RingFinder<ConnectTable>::ResultPointer>::const_iterator 
			it = blocks.begin(); it != blocks.end(); ++it)
		{
			Index nBonds = (*it)->countBonds();
			for (Index j = 0; j < nBonds; ++j) {
				const Bond& bondJ = (*it)->getBond(j);
				result[table.indexOf(bondJ)] = Bond::PARTIAL_BOND;
			}
		}
		return std::count(result, result+table.countBonds(), Bond::PARTIAL_BOND);
	}
	
}
